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Abstract 

The Melnikov criterion is used to examine a global homoclinic bifurcation and tran- 
sition to chaos in the case of a quarter car model excited kinematically by the 
road surface profile. By analyzing the potential an analytic expression is found for 
the homoclinic orbit. By introducing an harmonic excitation term and damping as 
perturbations, the critical Melnikov amplitude of the road surface profile is found, 
above which the system can vibrate chaotically. 

Keywords: Melnikov criterion, chaotic vibration, quarter-car, sky hook, magne- 
torheological dampers 



1 Introduction 



The problem of rough surface road profiles and its influence on vehicle un- 
wanted vibrations due to kinematic excitations is still a subject of research 
among automotive manufacturers and research groups, whose objective is to 
minimize their effects on the driver and passengers [1,2,3,4,5,6]. Past studies 
focused on the dynamics of a passive car suspension, the nonlinear character- 
istics of tyres and the effect of shimming in vehicle wheels [7,8,9]. Recently 
many new applications of active and semi active control procedures and special 
devices to minimize vehicle vibrations have been developed [10,11,12]. Conse- 
quently old mechanical quarter car models [8,9] have been re-examined in the 
context of active damper applications. Dampers based on magnetorheological 
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fluid with typical hysteretic characteristics have significant promise for effec- 
tive vibration damping in many applications [13,14,15]. New ideas in vehicle 
vibration damping, such as 'Sky hook' control [16] or control [15,17] have 
already been implemented and tested in several car and motorcycle applica- 
tions. Efforts have focused on studies of the excitation of the automobile by 
a road surface profile with harmful noise components [3,4]. However, due to 
various nonlinearities in the vehicle dynamics, chaotic behaviour may produce 
noise like responses [6,18,19]. 

In the present paper the model of Li et al. [6] is used, with the addition 
of a gravitational term that changes the equilibrium point and therefore the 
external potential. This paper uses the Melnikov theory [20,21,22] to estimate 
the critical amplitude of the road surface profile above which the system can 
vibrate chaoticaly. 

The gravitational term changes the topology of the heteroclinic orbit (in the 
case of a symmetric reversed 'Mexican hat' potential) into a homoclinic one 
(with a broken symmetry potential). Systems with Duffing characteristics hav- 
ing non-symmetric potentials and a linear repulsive force term are a wide class 
of mechanical systems and have been the subject of previous investigations in 
the context of the appearance of chaotic solutions [23,24,25,28]. This paper 
uses a similar approach for a quarter car model with a magnetorheological 
damper [13,14,15]. 



2 The quarter-car model 

The equation of motion of a single degree of freedom quarter-car model (Fig. 
1) is [6] 

d 2 x ( d \ 

m— + ki(x - x ) + mg + F h — (x - x ), x - x \ = 0, (1) 



where Fh is an additional nonlinear hysteretic suspension damping and stiff- 
ness force dependent on relative displacement and velocity, given by 



F h l^{x - x ),x - x j = k 2 (x - x ) 3 + ci-^ (x - x ) + c 2 (^(x-x ) 



(2) 



and 



xq = A sin(f2t). 



(3) 
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Defining a new variable for a relative displacement as 



y = X - X (4) 

we get 

+ u 2 y + B^ 3 + 5 2 ^f + B 3 (^J = - 5 + Afi 2 sin(m), (5) 

where a; 2 = ki/m, B\ = k 2 /m, B 2 = Ci/m, B 3 = c 2 /m. 
Following Li et al. [6] the system parameters are defined as 

m = 240 kg, h = 160000 N/m, k 2 = -300000 N/m 3 , (6) 
ci = -250 Ns/m, c 2 = 25 Ns 3 /m 3 . 



The corresponding dimensionless equation of motion can be written for a 
scaled time variable t — cut as: 

y + y + ky 3 + ay + [3y 3 = -g' + AVt' 2 sin(fiV), (7) 



where k = Bju 2 = a = B 2 /u = (3 = B 3 u = c 2 ^, g' = 

Q' = VL/uo. The overdots in Eq. (7) denote the corresponding derivative with 
respect to r ( = d/dr). 

Note that in our model, Eq. (5), we use both a complicated non-symmetric 
potential and also non-trivial damping of the Rayleigh type. Similar damping 
terms have been used before in the context of Melnikov theory [29,30]. Litak et 
al. [29] considered the Froude pendulum, with polynomial damping to model 
a dry friction phenomenon. Trueba et al. [30] performed systematic studies 
for basic nonlinear oscillators including those with combined damping. Here 
the motivation in using a complicated damping term is different, and arises 
from the use of magnetorheological dampers in vehicle suspensions [13,14,15]. 
The signs of the q coefficients (Eqs. 1-2) have changed compared to reference 
[6], in order to recover the usual Rayleigh term c 2 v 3 + civ, where v = dy/dt 
is the system velocity and c\ < 0, c 2 > (defined Eq. (6)). This term is able 
to drive the system into a stable limit cycle solution, being dissipative for a 
large enough velocity v (v > yjci/c^j and pumping energy for a small velocity 

(v < y/cijcz). 
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In Eq. (7), the nonlinear stiffness force has the potential 

V(y) = g'y + \y 2 + \y\ (8) 

Figure 2 shows this potential, and highlights the characteristic fixed points. 
Note the non-symmetry is caused by the gravitational term g'y, and that 
k < 0. 

In Figs. 3a-b we show the results of calculations in the interesting region of 
the main resonance for the system parameters given in Eq. (6) and a realistic 
amplitude of road profile excitation, namely A = 0.11 m. In this case the 
vehicle vibration amplitude, A OUT , plotted in Fig. 3a, has been determined 
numerically. For simplicity it has been defined as: 

A UT = (y m ax ~ 2/mm)/2, (9) 

where y max and y min are the maximum and minimum response of the vehicle 
model in the steady state. The resonance curve was calculated by tracking the 
solution for decreasing Q', and indicates that the main resonance occurs at 
Q' ~ 0.85. The response curve is inclined to the left, as expected for a nonlinear 
system with softening stiffness characteristic. Clearly, a jump between large 
and small vibration amplitudes exists at Q' ~ 0.8. Below this frequency we also 
observe a second, but much smaller, maximum of Aqut (at Q' ~ 0.75) which 
indicates that something interesting is occurring at this frequency. To examine 
this effect in more detail Fig. 3b shows a bifurcation diagram over the same 
range of excitation frequencies. Interestingly Q' m 0.75 is a point of dramatic 
change in the system behaviour. For any frequency below this point we see a 
black bounded region while above it there are singular points. One can easily 
see that the local change in Aqut is associated with a Hopf bifurcation. This 
transition is usually connected with a synchronization phenomenon between 
the system vibration frequency and the external excitation frequency. It is also 
connected with a slight change of the size of the attractor, reflected in the plot 
of A OUT (Q') (Fig. 3a). 

To show the changes in the dynamics caused by this resonance and Hopf 
bifurcation, Figs. 4a-c show the phase portraits and Poincare maps of the 
system for three chosen frequencies. One can easily identify Fig. 4a (Q' = 0.6) 
as a quasi-periodic solution with a limit cycle attractor. On the other hand 
the solutions presented in Fig. 4b (Q' = 0.8) and Fig. 4c (fl' = 1.1) show 
synchronized motion represented by singular points. The range of the vibration 
amplitudes is the highest for the last case examined (Q' = 1.1), which is 
consistent with Fig. 3a. 
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Figure 4d shows the hysteresis of the function Fh(y, y), defined by Eq. (2) and 
obtained during the same simulation sweeps as the phase portraits. The cor- 
responding hysteresis loops differ in size. Starting with T' plotted for Q' = 0.6 
then increasing strongly in size ('2' plotted for £1' = 0.8) and finally decreasing 
('3' plotted for f2' = f.l). Note this sequence differs from the changes in the 
vibration amplitude, where Q' = 1.1 has the largest amplitude, Aout, and 
this arises because of the combined effect of displacement and velocity. 



3 Melnikov Analysis 



Melnikov analysis starts with the renormalisation of the potential (Eq. (8), 
Fig. 2) [28]. If we let y = z + yo, where y is the fixed point given in Fig. 2, 
and V 1 (z) = V(y)-V(y ), then, 

V 1 (z) = jz 2 (z-z 1 )(z-z 2 ), (10) 



where Z\ = 1.298 and z 2 = 1.593. Fig. 5 shows this normalized potential. 
Notice that the left peak (the saddle point) of the potential Vi(z) occurs at 
z = < z 1 < z 2 and that Vi(0) = 0. 

Looking for a homoclinic orbit we introduce a small parameter e (formally 
ea = a, ef3 = (5 and eA — A). The equation of motion then has the following 
form, 

/ 3 1 \ 

z + eaz + e[3z 3 + k[z 3 - -(z x + z 2 )z 2 + -z x z 2 z = eAQ' 2 sin (fi'r). (11) 



Rewriting this second order differential equation as two first order differential 
equations yields, 



z = v, (12) 

3k k ~ 
v = -kz 3 + — (zi + z 2 )z 2 z x z 2 z + e(-av - (3v 3 + Ail' 2 sin (QV)). 



Note that the unperturbed equations (for e = 0) can be obtained from the 
gradients of the Hamiltonian H°(z,v), 

dH° dH° . . 

Z = ^ V = -^ (13) 
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where H° is defined as 



~2 + 1^ z ~ Zl >^ z ~ Z2 ' Z ■ 



(14) 



The homoclinic orbits are obtained from the unperturbed Hamiltonian, Eq. 

(14), as 



dz 



Zyf(z ~ Z X )[Z - Z 2 ) 



(15) 



which may be evaluated in the following form: 

2z]_z 2 - (zi + z 2 )z + 2^Jz 1 z 2 (z - zi)(z - z 2 ) 



T-T 



-ziz 2 k 



In 



, (16) 



where r is a time like constant of integration. 

Thus, the single homoclinic orbit is given by the inverse of the above expression 
and the corresponding velocity (z*(t — r ),v*(r — r )), as [28], 



4zi2 2 exp((T -r )J^f^ 



-( Zl - z 2 f - exp (2(r - r ) ^f^) + 2(zi + z 2 ) exp ((r - r )^f^) 



(17) 



■4 2l ^^exp ((r - r )^f^) ((^ - z 2 f - exp (2(* - t )^f^)) 



(-(*! - ^ 2 ) 2 - exp (2(r - r )^f^) + 2( Zl + z 2 ) exp ((r - r o ; 



-fcziZ2 



Now suppose that 



r = t i + r 02 , where r i = - In 



\[2{z 2 - zi) 



(18) 



t i has been fixed to guarantee the proper parity (under the time transforma- 
tion r — > — r), and hence 



z(—t) = z{t) and r) = — u(t). 



(19) 



r 2 is an arbitrary constant to be determined later in the minimization of the 
Melnikov integral M(t 02 ). The corresponding orbit (z*(r — t ),v*(t — r )) is 
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plotted on the phase plane in Fig. 6. 



The distance between perturbed stable and unstable manifolds and their possi- 
ble cross-sections may be examined by means of the Melnikov integral M(r 02 ), 
given by [21] 



M(r 02 ) = J h(z*(r - t 01 - t 02 ),v*(t - r i - r 02 )) 

— oo 

Ag(V (r - toi - Tea), v*(t - r i - r 02 ))dr, 



(20) 



where the wedge product for two dimensional vectors is defined as h A g = 
h\g 2 — h 2 gi. The corresponding vector h is the gradient of unperturbed Hamil- 
tonian (Eq. 13), 



H-z* 3 + -{z x + z 2 )z* 2 - ^ Zl z 2 z*),v* 



(21) 



while the vector g consists of the perturbation terms to the same Hamiltonian 
(Eq. (10)), 



-&v* - (3v* 3 + Att' 2 sin (fi'r), 



(22) 



Thus, shifting the time coordinate r — > r + r 02 under the integral (Eq. (20)), 
gives 



M(r 02 )= J v *(r-r 01 )((-av*(r-T 01 )-Pv* 3 (T-T 01 ) 

— oo 

+ AQ' 2 sm{Q'(T + T 02 )))dr. (23) 



Finally, a sufficient condition for a global homoclinic transition corresponding 
to a horseshoe type of stable and unstable manifold cross-section (for the 
excitation amplitude A > A c ), can be written as: 



V MM = and + 0. 



T02 



dr. 



(24) 



02 



From Eqs. (23) and (24) 



Ac ~wi 2 \n>y (25) 
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where 



h = 



(26) 



and 



I 2 (ft') = sup 

t 02 6R 



J v*(t - r i) sin (ft'(r + r 02 ))dr 

-00 

y v *(r — Tqi) sin (fi'r)dr 



(27) 



where sup means supremum for various r 02 , and is practically realized by 

cos(fi'r 02 ) = ±1. (28) 



Equation 27 has been obtained by using a trigonometric identity: sin('0 1 + 
■0 2 ) = sin(^i) cos (^2) + cos(V'i) sin(^ 2 ) where ^1 = Jl'r and ^ 2 = fi'r 02 . We 
left only the term cos('0 2 ) (Eq. 28) because of the odd parity of the velocity 
function v*(t — r i) under the integral. Of course to do such a simplification 
we needed defined parity of v* (Eq. 20) and a proper choice of a constant r i 
(Eq. 19). 

Note, the above integrals may be evaluated analytically [31] but here, for 
simplicity, they are calculated numerically [23,28]. Figure 7 shows A c as a 
function of ft' for a = —0.04 and (3 = 2.69 (see Eqs. (6) and (7)) given by 
the curve labelled '1' and f3 = 2.69/2 given by the curve labelled '2'. One can 
see the characteristic double sack-like shape, similar to the structure observed 
by Lenci and Rega [25]. This structure is governed by the oscillating term 
sin (Q't) in the denominator of the integral 7 2 (n') (Eq. (25)). 



4 Results of Simulations 



To illustrate the influence of a global homoclinic transition on the system 
dynamics, simulations were performed for interesting values of the system pa- 
rameters, using Eq. (7) for the model in Fig. 1. Knowing the critical value 
of the road profile amplitude A c (Fig. 7) and looking at a typical homoclinic 
bifurcation [22] the effect on the resonance curves were examined first. Fig- 
ure 8 shows the sequence of resonance curves for A = 0.11, 0.16, 0.21, 0.26, 
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0.31 and 0.36m respectively. Apart from a typical shift of the maximum re- 
sponse to the right all of these curves are very similar, up to A = 0.31m. 
For A = 0.41m the synchronized solution is not stable in the region of reso- 
nance. The other difference in the system behaviour occurs to the left side of 
the resonance peak where multiple solutions of the nonlinear system appear 
(in this case resonant and non-resonant solutions) in the region of the reso- 
nance. One can observe that starting from A = 0.26m the curves in Fig. 8 
show a discontinuity signaling jumps between the resonant and non-resonant 
vibration amplitude A OUT . Note in all cases a series of simulations were per- 
formed to calculate the system response, with Q' decreasing as in Figs. 3a-b. 
For most of curves the same initial conditions were used for large Q', namely 
[xj n , Uj n ] = [0.15,0.1]. However if the system escaped from the potential well 
initial conditions of [x in ,v in ] = [—0.15,0.1] and [0,0.1] were used to avoid this 
effect. For A = 0.36m this was not possible in the vicinity of the resonance 
peak where the system escaped from the potential well for any initial con- 
ditions. Moreover just before this escape (for A = 0.36m m A c ) we observe 
a further increase in the vibration amplitude Aout- Examining the related 
bifurcations diagrams we identified period doubling phenomenon occurring in 
this region which can be classified as a precursor of chaotic vibrations. Indeed 
alternative criteria to the Melnikov approach (Eq. 24) are based on the period 
doubling cascade [26,27]. For larger amplitudes the unstable vibration region 
where escape from the potential is possible increases. On the other hand, at 
A = A c the border between the basins of attraction belonging to different 
solutions disappears. To avoid these difficulties for further analysis the syn- 
chronized solution for Q' = 0.8 at A = 0.31m (Fig. 8) was used, and then the 
excitation amplitude was increased slightly to A = 0.41m, crossing the criti- 
cal amplitude of A c = 0.35m. Figures 9a-b show the phase portraits (by lines) 
and Poincare maps (by points) for these two cases. Figure 9a shows a synchro- 
nized motion while Fig. 9b corresponds to a chaotic attractor. The dominant 
Lyapunov exponents calculated for these responses are Ai = —0.1625 (Fig. 
9a) and Ai = 0.03540 (Fig. 9b). Note, the chaotic attractor is very similar to 
that studied by Thompson [32] where the harmonic potential has been sup- 
plemented by a nonlinear term with displacement to the power 3 {z 3 ). Figure 
9 also shows the time histories for the two cases: A = 0.31m (Fig. 9c) and 
A = 0.41m (Fig. 9d). In this figure the difference between the periodic and 
chaotic responses is clear. 



5 Summary and Conclusions 

We have studied the vibrations of a quarter-car model with a softening stiff- 
ness of the Duffing type, focusing on the potential for chaotic behaviour. The 
model and parameters used were taken from the paper by Li et ai [6], with 
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the addition of the gravity force. The addition of this gravity force breaks 
the symmetry of the potential, so that V{— x) 7^ V(x). The hysteretic na- 
ture of the damper caused a range of interesting system behaviour, such as 
quasi-periodic, synchronized and chaotic motion. This had a substantial effect 
on the heteroclinic orbits, which transformed into homoclinic orbits. We ex- 
amined the global homoclinic bifurcations that appear as instabilities at the 
boundaries of the basins of attraction, and the cross-sections of stable and 
unstable manifolds, by the perturbation approach using Melnikov theory. A 
critical amplitude was found for which the system can exhibit chaotic vibra- 
tions. The analytic results have been confirmed by numerical simulations. In 
particular, the chaotic strange attractor was found for an excitation amplitude 
A at the critical value, A c , and a period doubling precursor for A = 0.36m. 
The transition to chaos appears to be present for A c = 0.36m but could be low- 
ered significantly for a smaller damping coefficient c 2 (Fig. 7). Fortunately this 
region is beyond the usual amplitude of road profile excitation. The chaotic 
solution appears just before the escape from the potential well, which is similar 
to the system with a non-symmetric potential described by Thompson [32]. 
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Fig. 1. The quarter-car model subjected to kinematic excitation with nonlinear 
damping and stiffness. 




— 0.08 



0.04 



0.00 




Fig. 2. External potential V(y) (Eq. (8)) for given system parameters (Eq. (6)). V(y) 
is scaled in Nm while y is in m. The fixed points are yo = —0.7228m, y\ = —0.0147m, 
2/2 = 0.7375m. 
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Fig. 3. Vibration amplitude Aout = (Umax — Umin)/^ (Fig. 3a) and bifurcation 
diagram (Fig. 3b). The amplitude of a road profile has been taken as A = 0.11m. The 
arrows indicates reduces in the simulations. For each new smaller the initial 
conditions [yi m Vi n ] were taken as the final position and velocity for the previous Q' . 
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Fig. 4. Phase diagrams (velocity v = y versus displacement y plotted by lines) and 
corresponding Poincare sections (plotted by points) for A = 0.11m and different 
W: Q' = 0.6 (Fig. 4a), Q' = 0.8 (Fig. 4b) and Q,' = 1.1 (4c). The corresponding 
hysteresis curves are shown in Fig. 4d, where '1', '2' and '3' represent O' = 0.6, 
0.8 and 1.1, respectively, vuj is scaled in m/s, y in m, while the renormalized is 
presented in dimensionless units Fh = ay + (3y 3 + ky 3 (see Eq. (7)). 
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0.04 




Fig. 5. Renormalized external potential V\{z) = -^z (z— z\){z — z%) for given system 
parameters (in our case: k = —1.875 N/m 3 , z\ ~ 1.298m and Z2 ~ 1.593m). z is 
expressed in m while V(z) is in Nm. 
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Fig. 6. A homoclinic orbit for the given potential (Eq. (9) and Fig. 5). z is expressed 
in m while vuj is given in m/s. 



16 



'2' 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 



Fig. 7. Critical amplitude A c versus vibration frequency for two different damping 
choice. The damping coefficients C2 = 25 Ns 3 /m 3 (for other system parameters see 
Eq. (6)) , and c 2 = 12.5 Ns 3 /m 3 , for '1' and '2' curves respectively. The points 'reg' 
and 'cha' denote regular and chaotic solutions (to be examined later in Figs. 9a and 
b). 
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Fig. 8. Sequence of vehicle vibration amplitudes Aout = (Umax — Umin) /2 versus 
frequency Q' for road profile amplitudes of A = 0.11, 0.16, 0.21, 0.26, 0.31 and 
0.36m (from the lower to upper curves, respectively). 
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Fig. 9. Phase portraits and Poincare maps below and above the critical amplitude 
at = 0.8 for regular (Fig. 9a, A = 0.31m and initial conditions [xi n ,Vi n ] = 
[0.15,0.1]) and chaotic (Fig. 9b, A = 0.41m and the initial conditions [xi n ,Vi n ] = 
[—0.15,0.1]) solutions. The dominant Lyapunov exponents are Ai = —0.1625 and 
Ai = 0.0354, respectively. Figs. 9c and 9d give the corresponding time responses. 
Here vlo is given in m/s while y is in m. 
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